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Abstract 

We extend our previous numerical simulation of accretion disks with shock 
waves when cooling effects are also included. We consider bremsstrahlung and 
other power law processes: A oc T a p 2 to mimic cooling in our simulation. We 
employ Smoothed Particle Hydrodynamics technique as in the past. We observe 
that for a given angular momentum of the flow, the shock wave undergoes a steady, 
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radial oscillation with the period is roughly equal to the cooling time. Oscillations 
seem to take place when the disk and cooling parameters (i.e., accretion rate, 
cooling process) are such that the infall time from shock is of the same order 
as the post-shock cooling time. The amplitude of oscillation could be up to ten 
percent of the distance of the shock wave from the black hole when the black 
hole is accreting. When the accretion is impossible due to the centrifugal barrier, 
the amplitude variation could be much larger. Due to the oscillation, the energy 
output from the disk is also seen to vary quasi-periodically. We believe that 
these oscillations might be responsible for the quasi periodic oscillation (QPO) 
behaviors seen in several black hole candidates, in neutron star systems as well 
as dwarf novae outbursts such as SS Cygni and VW Hyi. 
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1. INTRODUCTION 



It is becoming increasingly clear that if in active galaxies accretion takes place 
in the form of accretion disks, then these disks may not be similar to the standard 
'Shakura-Sunyaev' type thin accretion disks models found most suitable for binary 
systems. The presence of virtually simultaneous variabilities in optical and UV in 
several objects such as NGC 5548 (Clavel et al. 1990, 1992 and Peterson et al. 1991), 
NGC 4151 (Perola et al. 1986) and temporal variation of the line emissions in broad line 
radio galaxies, such as ARP 102B, 3C390.3 (Miller & Peterson, 1990; Veilleux & Zheng, 
1991; Zheng, Veilleux & Grandi, 1991; Chakrabarti & Wiita, 1994) are examples of the 
fact that the some disks are of completely different type and a significant contribution 
to the emitted spectra could be made by the shock waves in the disks. These shock 
waves could either be pressure and angular momentum supported (Chakrabarti, 1989 
hereafter C89; Chakrabarti, 1990 hereafter C90) in the case of inflows which include 
low but significant angular momentum, or only pressure supported (Chang & Ostriker, 
1985; Kazanas & Ellison, 1986; Babul, Ostriker & Meszaros, 1993) when the flow is 
spherical. Whether or not accretion disks may have a shock or simply accrete onto the 
black hole through inner sonic point could be determined from the viscosity, accretion 
rate and angular momentum of the flow. This unified classification of the global solution 
is discussed in C89. 

In our previous work (Chakrabarti & Molteni 1993; hereafter CM93), we pre- 
sented hydrodynamic simulation of the standing shock waves in thin accretion disks. We 
showed in particular, that standing shock waves are formed in inviscid accretion disks 
of low angular momentum A < A m 6, where X m b is the marginally bound angular mo- 
mentum. Subsequently, Molteni, Lanzafame & Chakrabarti (1994; hereafter MLC94), 
studied numerical simulation of shock waves in quasi-spherical accretion disks around a 
Schwarzschild black hole and verified that the results agree very well with the vertically 

3 



averaged solution (C89). Sponholz and Molteni (1994) studied shock formation around 
a Kerr black hole and find differing shock locations in co- and contra-rotating flows. 
When angular momentum supply at the outer edge of the disk is large, a significant 
viscosity must be present in the disk in order that accretion may take place. In an 
important simulation, Chakrabarti & Molteni (1995, hereafter CM95) pointed out that 
when viscosity is low enough (in the language of a v i s parameter of Shakura-Sunyaev 
1973; a v i s < a vc ) } a weaker shock may still form farther away from the black hole. For 
a higher viscosity (a v i s > a vc ) } the disk is shock-free and is similar to a standard Kep- 
lerian disk at larger distance while passing through an inner sonic point. These results 
are consistent with the observation of the most general solution of the viscous transonic 
disks (C90) where the critical viscosity parameter (a vc ) for isothermal viscous disk is 
computed. Recent works of 'newly discovered advection dominated' flows (Narayan & 
Yi, 1994, Narayan & Yi, 1995), which are the generalization of shock-free global so- 
lutions of C90, (for cooling processes not satisfying isothermality condition) miss the 
shock solutions completely because of the restrictive assumption of self-similarity and 
the choice of the Newtonian potential. Thus the observations derived from our more 
general consideration will not be applicable to such solutions. 

In the present paper, we deviate from our previous studies of simulation in two 
ways: (a) we use an adiabatic index 7 = 5/3 (instead of 7 = 4/3) to mimic a gas with 
a low radiation pressure and (b) we include a a power-law cooling process A oc p 2 T a in 
the disk. We extensively study the case of a = 0.5 (i.e., bremsstrahlung) and discover 
a completely new physical effect not discussed so far in the literature in the context of 
transonic accretion disks around a black hole. We show that the shock wave undergoes 
a steady quasi-periodic oscillation in radial direction with an amplitude roughly about 
five to ten percent of the steady state shock distance from the black hole and the 
time period similar to the cooling time. The oscillation causes the luminosity of the 
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disk to vary quasi-periodically as well. We find a resonance behavior in the oscillation 
characteristics in the accreting disk: if the disk parameters (such as the accretion rate 
and angular momentum) and the cooling parameters (such as a and efficiency of cooling 
relating to heating) are such that the infall time from the shock location is comparable 
to the radiative cooling time, then the shock oscillates, otherwise only the steady state 
shock remains. 

One of the motivations of our study stems from the results of the Langer, Chan- 
mugam & Shaviv (1981, 1982) who showed that shocks in accretion columns of white 
dwarfs and neutron stars exhibit oscillational instability in presence of a power law cool- 
ing. Subsequent work of linear stability analysis (Chevalier & Imamura 1982) shows 
that such an oscillation is expected. These results were further generalized by several 
workers by including effects of cyclotron cooling (Chanmugam, Langer, & Shaviv, 1985) 
in accretion onto magnetized white dwarfs who found that the cooling at a field strength 
B ~ 30 mG is sufficient to stabilize the oscillation. Wolff, Gardner & Wood (1989) 
studies one and two temperature models including bremsstrahlung, Compton cooling 
and cyclotron cooling and Wu, Chanmugam & Shaviv (1992) studies oscillations in con- 
tinuously excited accretion flows with a boundary condition relevant for white dwarfs. 
Other works discussed analogies between the instabilities in shocks due to cooling and 
exothermic reactions (e.g. Falle, 1995 and references therein). Our results are new in 
the following respects: Contrary to white dwarf or neutron star, a black hole does not 
have a 'hard' surface. However, in presence of weaker dissipation, the centrifugal barrier 
experienced by the flow can be strong enough to act as a 'hard' surface where eventually 
a standing shock wave is formed. Thus, the question arises: would these shock waves 
oscillate in the manner it does on white dwarf surface? We find that they do and we 
believe that they may be responsible for the QPO behaviors observed in many black 
hole candidates. Secondly, we found that the oscillation is actually of 'resonance' type, 
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namely, it occurs only the the cooling time scales roughly match with the infall time 
scale. If the cooling time is very short the flow rapidly adjusts to the cooler disk, but if 
the cooling time is very large, the thermal energy is advected along with the flow. As 
in earlier three-dimensional simulation (Molteni, Lanzafame & Chakrabarti, 1994), we 
find that the post-shock region forms a corona which could be the site of hard X-ray 
and 7-ray production. 

The general theory of the study of the shock waves in accretion disks are presented 
in C89, C90 and CM93 and we shall not repeat them here. In the next Section, we 
shall present the equations used in our numerical simulations. In §3, we present the 
way bremsstrahlung is implemented in our axisymmetric SPH code. In §4, we present 
numerical solutions and compare with the analytical results. Finally, in §5, we make 
concluding remarks. 

2. MODEL EQUATIONS 

We assume a rotating, axisymmetric, accretion flow around a black hole. We 
take Newtonian model for the non-rotating central compact object as given in terms 
of the Paczyhski & Wiita (1980) potential which is found to be accurate enough for 
astrophysical purposes. The internal energy e of the disk is defined as P = pe(-f — 1). 
P and p are the isotropic pressure and the matter density respectively, 7 = 5/3 is the 
adiabatic index. Only single temperature model is considered here so that electrons and 
protons have the same temperature (two temperature disk solutions will be presented 
elsewhere). We presently assume that the flow is inviscid, i.e., the specific angular 
momentum A is constant everywhere. We use the density of the disk at the outer 
edge p re j = p 0} the velocity of light v re j = c and the Schwarzschild radius x re j = 
R g = 2GM/c 2 of the black hole mass M, as the reference density, velocity and distance 
respectively. Thus, for example, the unit of time is t re j = x re j/v re j = R g /c } the unit 
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of angular momentum is X re j = x re jv re j = cR g} the unit of mass is M re j = p r ef%ref = 
p re jR^ } and the unit of mass accretion rate is M re j = p re jx^ e j /t re j = p re jR 2 c. 

First, we consider one-dimensional equations. Our choice is restricted by the fol- 
lowing constraints. Whereas it is easy to obtain analytical solutions for disks of constant 
thickness (e.g., CM93) or disks in vertical equilibrium (e.g., C89), it is very difficult to 
obtain solutions for full three-dimensional flows. An axisymmetric numerical code can 
easily test analytical solution of one-dimensional axisymmetric flow of constant thick- 
ness and it can also perform a fully three-dimensional simulation without taking resort 
to the vertical equilibrium assumption. Thus, we solve equations of constant thickness 
(rather than vertically averaged model) first to verify that the time dependent code is 
satisfactory. Subsequently, we solve fully three-dimensional equations as presented in 
§3. The time-dependent equations which we solve are: 

(a) The energy equation: 

de de P dv 

dt dx p dx 1/,2/? 

(b) The radial momentum equation: 

dv dv 1 dP , ., N 
at ox p op 



(la) 



(16) 



(c) The continuity equation: 

dp 1 d . . . 

i + ^ m) = • 

Here, x is the radial coordinate, ip(x) is the effective potential energy ip(x) = g(x) + 
A 2 /2x 2 with g(x) as the radial force potential, which, in the pseudo-Newtonian model, 
takes the form: g(x) = —l/2(x — e is the internal energy density e = p^yy^, v is 

the radial component of velocity and a is the adiabatic sound speed: a 2 = ^ Primes 
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denote a derivative with respect to x. ( a is the non-dimensional bremsstrahlung loss 
coefficient, 

f T l J 2 



an> 



d 



Jp re fX re fl re j 
Cl/2 = — 



Tref = Pj ^f 1 (3) 



where, p = 0.5 and j = 1.4 X 10 2 c.g.s. unit for ionized hydrogen (Allen 1973), 



m r 



is the mass of the proton, and k is the Boltzmann constant. The subscript 1/2 of ( 
reflects that we are using the cooling law A = Ci/2P 2 T a with the constant d/ 2 exactly 
same as in bremsstrahlung case (a = 0.5). 

In the steady state, the above equations could be solved quite easily in one di- 
mension by assuming d/dt = and re-arranging them first in the form: 

« = dx x 3^ ( 4a ) 

ax a z jv — v 

which along with the following equations, 

dv 2a da a 2 dp dib , , . 

ax 7 ax ax ax 



an> 



d 



1 dv 1 dp 1 , . 

-T + -r + - = °- ( 4c ) 

v ax p ax x 

Detailed solution properties of the sonic points and shocks in presence of cooling effects 
are beyond the scope of this paper and will be discussed elsewhere. 

A complete solution which includes a shock wave 'glues' together two integral 
curves passing through the outer and inner sonic points respectively at the shock loca- 
tion x s determined by the Rankine-Hugoniot relation at the shock (Chakrabarti, 1990). 
Assuming the shock to be infinitesimally thin and non-dissipative at the shock itself, 
these relations are of (a) the conservation of the energy flux, (b) the conservation of the 
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momentum flux and (c) the conservation of the mass flux, respectively across the shock 
wave. When these relations are combined together, one obtains the so-called 'Mach 
number relation' as: 

_J£L±i_ = K + 3 = Constant (5) 

where, — and + signs are used to denote quantities in the pre-shock and post-shock 
regions. This relation is used to obtain the location of the shock wave in a steady state. 
Before we present the analytical and numerical solutions, we present a few sentences 
about the way the cooling effect is implemented in smoothed particle hydrodynamics 
simulations (SPH). 

3. SPH: IMPLEMENTATION OF THE COOLING LAW 



To find non-steady-state solutions of the classical one-dimensional hydrodynam- 
ical equations l(a-c) of a compressible inviscid fluid we apply the Smoothed Particle 
Hydrodynamics (SPH) method, specially adopted for axial symmetric configurations as 
in our earlier works. The procedure presented here is valid for thick (r, z) accretion as 
well and is described in detail in Molteni & Sponholz (1994) as well as in CM93. While 
applying to the thin disks (la-c), one has to assume a constant vertical thickness and 
inject matter only on the equatorial plane. 

The basic point for the cylindric geometry approach is to assume the interpolating 
Kernel W which is a function of cylindrical radial coordinate r = r(x,z) and k-th 
particle of mass m k as, 

m k = 2irp k x k Ar k (6) 
Any smooth function A(fi) at r 4 - is defined as, 

A(f t ) = ( A(r)W(r - f t ; Ky^-dr w ]T -^—A{r k )W{r k - f t ;h), (7) 
J Znpx ^ ZTTX k p k 
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h being the particle size. Thus, for example, for the density at each point, we have 
the simple expression that identically satisfy the continuity equation in the cylindrical 
coordinate: 

p{?i)*Y,— W ( ? k-?i\h) (8) 

k Xk 

The equations of motion to be solved consist of the radial momentum equation: 

Dv x \ _ / El ,E± „ \ dW tk P _ 1 Xj 

DtJr \ x k \pl + pl + Ui V d Xi + xf 2(r,- - l) 2 r,- ' W 

the vertical component of the momentum equation. 

DVz \ _ ^ ( El + E i i n ^ dW ik 1 Z i (Q h \ 

DtJr \ x k \pl + pl + Ui V dzi 2(r,- - l) 2 r,- " l " j 



The energy equation describes the behavior of the thermal energy per unit mass 
and contains the cooling-term A; = Ci/2Pi(£i) a '- 

\DtJi %k X k \Pi Pk J Pi 

To mimic numerically the kinematic dissipation, the artificial viscosities: 



Pij 



~ ^i^xi ^j^xj zi ^zi){^i %j) 

1113 = *m + 1 2 ) w+V 2 ' 

l\ = (x t - Xjf + (z t - zjf, r] = O.lh 2 
with the abbreviations for density p^ and sound speed a 8J are used (Monaghan 1992): 

Pi + Pj t _ a, + a,j 
Pa = — g — aij = — 2 — 

r] = 0.1h 2 , 

a and (3 control the amount of the artificial viscosity, necessary to reduce oscillations 
in shock transitions. 
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Thus, the implementation of the SPH in cylindrical code is similar in every respect 
to the original SPH code in Cartesian coordinate, except that the mass of a particle 
appears divided by its axial distance and the relative velocity — V{ between &-th and 
z-th particles must be replaced by a (x^Vk — X{Vi)/xi term (Chakrabarti & Molteni, 
1993; Molteni k Sponholz, 1994). 

4. RESULTS AND PHYSICAL INTERPRETATIONS 

Below, we first present analytical and numerical solutions in thin accretion disks 
with constant thickness. Subsequently, we present more realistic solutions for three- 
dimensional accretion disks. In the next sub-Section, we provide physical interpretation 
of these solutions. 

4.1 Results in Thin Accretion Disks 

While computing the results we have fixed our attention keeping an application 
to the active galaxies in mind. We have chosen the mass of the black hole to be 
M = 1O 8 M . We choose the specific angular momentum (in units of 2GM/c) to be 
A = 1.9 which lies between marginally stable and marginally bound values. We inject 
matter at the outer edge of the disk x = 50R g with v ~ 0.14 and a ~ 0.0024. Matter 
reaching the inner boundary (chosen here to be at x = 1.5) are removed. 

To obtain an analytical solution, we choose a typical set of parameters (such as, 
the density /?, the velocity v, the sound speed a and the specific angular momentum 
A) at the boundary and solve for p(x), v(x) and a(x) using equations 4(a-c). We have 
obtained steady state solutions using fourth order Runge-Kutta equations. 

Figure 1 compares the shock locations as a function of time for different value 
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of a. In all the cases, the density at the outer edge as 1.85 X 10~ 15 gm cm -3 . At 
this density, with a uniform vertical thickness of one Schwarzschild radius, the mass 
accretion rate is given by M = O.14M yr -1 which is about seventy percent of of the 
Eddington rate for a black hole of mass 1O 8 M . The solid curve is for a = 0.5, and the 
long dashed and short dashed curves are for a = 0.4 (largest amplitude) and a = 0.75 
(smallest amplitude) respectively. A large number of important conclusions could be 
drawn: (a) for bremsstrahlung cooling (a = 0.5), the shock oscillates in a very stable 
manner. The shock location x s varies by about five percent of the mean location. The 
period of oscillation is Ti/ 2 ~ 1500 unit, which is equal to ~ 1.5 X 10 6 s or 17. 4d. (b) For 
a = 0.75, the oscillation of the shock damps out completely, and the mean location is 
higher than the average x s for a = 0.5. In both (a) and (b), matter is found to accrete 
through the inner sonic point, (c) For a = 0.4, the amplitude is very high. In this case, 
no accretion takes place. The flow in the post-shock region cools down so rapidly that 
the pressure is not sufficient to push matter across the centrifugal barrier. Presence 
of viscosity is expected to reduce this barrier as well as the amplitude of oscillation. 
The time period of oscillation is about T p ~ 1978 unit, i.e., 22. 9d for a 10 8 M e black 
hole. In all these simulations the numerical length-scale (interpolation-length h) size is 
chosen to be 0.25 which is sufficiently small compared to the amplitude of oscillation. 
Therefore we believe that these oscillations are physical. The number of particles with 
the region of integration was about 200 in our simulation. 

In order to show that these oscillations are not artifacts of our simulation, we 
have tested the result in numerous ways, particularly by varying the particle size h and 
viscosity parameters. We find no difference of any significance. More importantly, we 
tested these solutions against analytical steady state solution using same method as in 
C89 and C90. In Fig. 2 we present one such test result where both the time dependent 
(Eqs. la-c) and time independent (Eqs. 4a-c) solutions are shown. We have chosen here 
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a = 0.5. The solid curves (supersonic and subsonic branches) represent the steady state 
result, while the long-dashed, short-dashed and dotted curves are the time dependent 
solutions at different phases of oscillation. The time dependent solutions, particularly 
in the subsonic branch, are seen to oscillate around the steady state solution. For 
comparison, we present the dash-dotted curve which is the steady state solution with 
no cooling effects present. In the pre-shock region, temperature and densities are low 
and cooling effects are negligible. In the post-shock region, cooling effects are very 
strong. In Fig. 3(a-b) we show how the (a) temperature (in units of 10 9 K) and (b) 
total bremsstrahlung loss (in arbitrary units) vary with distance from the black hole, 
as well as at various phases of oscillation. It is important to note that the whereas the 
temperature varies significantly throughout the post-shock region, maximum variation 
of the loss takes place much closer to the black hole at X ~ 5, because of higher 
densities (Of course, significant variation occurs at the mean location of the shock, due 
to the movement of the shock). When the shock is at its outer most position (short 
dashed curve) the bremsstrahlung loss in minimum. However, the net loss is still a 
strongly varying function of time. Figure 4 shows the variation of net loss with time. 
The solid curve is for a = 0.5 and the dashed curve is for a = 0.75. The variation of 
bremsstrahlung luminosity is about 15 percent. We do not find any steady accreting 
solution for the above accretion rate (i.e. M = O.GMEdd) when a = 0.4. The time 
dependent solution shows larger variations in the amplitude of oscillation (Fig. 1). We 
shall comment on the implications of this results in next Sections. 

4.3 Simulations of Thick Accretion Disks 

So far, we have presented the behavior of shock waves in a thin accretion disk. 
Presently, we show behaviors of axisymmetric three dimension axisymmetric flows 
whose SPH implementation is given in Eqs. 9(a-c). In these simulations we have 
assumed the disk to be symmetric with respect to the equatorial plane. Solutions with- 
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out the latter boundary conditions show very interesting behavior and we shall discuss 
about them elsewhere (Sponholz, Molteni & Chakrabarti 1995, in preparation). In Fig. 
5(a-d), we show results of three dimensional simulations. The parameters at the outer 
edge (X = 50) are chosen to be: density p = 4 X 10~ 14 gm c -1 , infall velocity u = 0.126 
and the sound velocity a = 0.04. This corresponds to an accretion rate of M = 1O.8M 
yr -1 , or, about M ~ 47MEdd- The specific angular momentum chosen is A = 1.75. 
With this low angular momentum, matter is still optically thin except extremely close 
to the black hole which is much inside the shock region and are not expected to influence 
our solution. The number of particles in the simulation is about 2000. In the upper 
half of each panel we show velocity vectors (whose sizes are proportional to the actual 
velocities), the thick curves are the the contours of constant Mach number M = 1. The 
thin lines presents on both (sub- and supersonic) sides represents contours of constant 
Mach numbers log(M) = —1.0, —0.5, and 0.5, 1.0, 1.5, ...3., whereas the dashed contours 
mark the subsonic matter. In the lower half of each panel we show contours of constant 
density in units of the density at the outer edge. The circle around the origin is of three 
Schwarzschild radius, the size of the marginally stable orbit. In Fig. 5a, we produce 
the solutions without bremsstrahlung. A very weak shock is formed very close to the 
black hole X ~ 5 and a strong wind is also produced from the high entropy, post-shock 
flows very similar to what is observed in MLC94. The solution roughly remains the 
same with time. 

When bremsstrahlung loss is included, the situation is changed dramatically. The 
shock is formed at around X = 10 and a strong corona with outflowing wind is also 
produced. As the shock oscillates, the size of the corona also changes. The net lumi- 
nosity due to bremsstrahlung effect is also varied by 15 to 20 percent. The time period 
of oscillation is about 2200. Figure 6, shows the oscillation of the emitted radiation 
from the three dimensional simulations. The oscillation at different region of the disk 
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differ slightly by frequency, which causes some smoothening effects, i.e., some oscillating 
region followed by some region of cancelation. 



4.3 Interpretation of the Results 

To understand the behavior of the shock waves qualitatively, we first note that 
the shocks discussed here are very strong, i.e., M_/M + >> 1. In this case, it is easily 
shown that, 

P+ = \p-vl. (10a) 

The pre-shock flow has a very low angular momentum and is almost freely falling. Thus, 
as the shock, 

X s 



and, 

1 , , 11 

x 

The density of the post-shock flow is 



v + (x s ) = -v_(x s ) = - — . (10c) 



4M , „ 

p+ ~ a p- = — ( 10o O 

X s 

and from Eqn. (10a), the sound speed at the post-shock position is obtained as, 

a + (x.s) = j(f ) 1/2 - (10e) 

x s 

The post-shock Mach number of the flow is, 

M+ = — = 0.2 1/2 ~ 0.45. (10/) 



at the shock. The location of the shock is determined by employing the Mach number 
relation (Eqn. 5) by noting that the shock invariant constant for a strong shock M_ — > 
oo is given by 

C = l/ 7 2 = 0.36. (lOflf) 
15 



In all our solutions we do see that the post-shock Mach number is ~ 0.45 (e.g., Fig. 2) 
and at the shock location C ~ 0.36 as described above. 

As the shock is perturbed and propagates outwards, it heats the post-shock flow 
to a higher temperature since the relative velocity (which is thermalized in the post- 
shock region) between the shock and the incoming flow is higher (Fig. 3a). As a result, 
the cooling rate goes up and the outward motion of the shock stops when the flow is 
sufficiently cool. The post-shock pressure drops and the shock returns towards the black 
hole. The relative motion between the shock and the post-shock flow decreases and the 
temperature in the post-shock region drops. The lower pressure in the post-shock region 
is unable to balance the pre-shock ram pressure and the shock collapse continues till the 
centrifugal barrier is sufficient to hold the flow. The shock then overshoots this region 
and bounces from higher pressure region and the cycle continues. In order to have a 
oscillatory behavior, the post-shock region must be able to cool in a time-scale t coo \ 
comparable to the advection time-scale t a dv If t coo i << t a d v , the flow quickly adjusts 
to the cooler state and if t coo \ >> t a d v the flow advects all the thermal energy. Only 
when these time scales are comparable, the flow tries to constantly adjust itself, and 
the oscillation is seen. Thus, the time period of oscillation t oscn is roughly comparable 
to the cooling time t coo \ ~ t a dv In above simulation, we have chosen the accretion rate 
in a manner that these time scales are comparable when a = 0.5. For this particular 
accretion accretion rate, t coo \ << t a d v , when a = 0.75 and the flow adjusts immediately 
to the new solution. As a result there is no oscillation. For a = 0.4, for the same 
accretion rate, no transonic accreting solution is possible. Here, the cooling time is 
large and the shock has to move outwards by a large amplitude (and for a longer time) 
so as to satisfy t coo \ ~ t oscn . The amplitude of oscillation is proportional to the thermal 
energy content of the flow. 
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To compare the time scales in our problem, we compute 




and the cooling time due to power-law cooling 




for various a. For a = 0.5, t a d v = 7.5 X 10 5 s and t coo \ = 1.35 X 10 6 s. The time period 
of oscillation (Fig. 1) is found to be t oscn = 1.5 X 10 6 s. Thus, t coo \ t oscn and t a d v 
is also similar within a factor of 2. On the other hand, for a = 0.75, we compute, 
tadv = 1-03 X 10 6 s and t coo \ = 3 X 3s. Because of the disparity in time scales, no 
oscillation is seen. For a = 0.4, cooling is very high and no transonic solution exists 
with the accretion rate chosen. The infall time-scale (t a( [ v ~ Ax^J 2 ) from x s ~ 42 to 
x s ~ 15 is given by t a d v = 865units or t a d v ~ 8.6 X 10 5 M/10 8 M e s roughly the same 
time in which half the oscillation is completed (Fig. 1). Thus, the oscillation is purely 
driven by the centrifugal barrier and the cooling effects: Centrifugal barrier determines 
the minimum shock location, whereas the amplitude of oscillation is determined by the 
time scales of various cooling processes. 



5. CONCLUDING REMARKS 



In this paper, we have discovered an important new physical effect in the context 
of the black hole accretion physics. We showed that if the disk parameters are right 
enough so that the infall timescale of accretion matter from the shock, roughly agrees 
with the cooling timescale in the post-shock region, then the shock exhibits a radial 
oscillation around the mean shock location with the time period as that of the cooling 
time scale. We have provided examples of such oscillation as well as the damping of the 
oscillation when one is farther away from 'resonance' behavior. We used a power law 
cooling process, such as bremsstrahlung, to illustrate our claim. For a given accretion 
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rate, there appears to be a cut-off exponent a c such that the oscillation is present when 
a < alpha C} however, it is possible that no 'universal' a c exists. If the shock forms for 
a given disk parameter, there could be always some a for which the oscillation of shock 
may be seen. 

Similar phenomenon has been reported in accretion processes onto compact stars 
which have 'hard' boundaries such as white dwarfs (Langer, Chanmugam & Shaviv, 
1981). It is found that a c 0.6 (Imamura, Wolff & Durisen, 1983). However, in these 
systems the location of the shock may form only very close to the star surface (unlike in 
the case of a black hole accretion where one can adjust the angular momentum to obtain 
the shock location far away as well), and therefore there may be a fixed a c as there is 
little room to vary the infall time scale. Secondly, we believe that the oscillation should 
take place for any cooling processes (not only with power law cooling process) provided 
the time scales agree. In other words, if an oscillation is seen for a given cooling process, 
it may disappear in presence of some other cooling process if the flow are not changed. 
This may be the reason for the disappearance of oscillation of Chanmugam, Langer & 
Shaviv (1985) when synchrotron cooling was added. 

The time period of oscillation is found to be of the order several hundred light 
crossing time of the black hole and therefore could range from milliseconds to months 
depending upon the mass of the black hole. The amplitude of the emitted radiation, 
mostly in regions extreme ultraviolet to 7 rays, should be modulated by as much as 
ten percent for accretion solution. Because the cool and hot matter co-exist side by 
side close to a shock, correlated variability is expected in these systems, which may 
have been observed in some cases as mentioned in the Introduction. The accretion 
rates needed are reasonable and we expect such phenomenon to be seen in many more 
systems, particularly in active galaxies where low angular momentum accretion may 
be taking place. The oscillation frequency only indirectly depends on viscosity since 
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the shock location is affected by viscosity (C90, CM95). In several galactic black hole 
candidates, such as GS339-4 and GS1124-68, QPO behavior is seen (Dotani, 1992). 
These oscillations are of 3 — 8Hz and the amplitude modulation could be as high as 
5 — 10 percent depending upon observed energy band. We know of no other physical 
explanation which can satisfactorily explain so large amplitude variation. Within the 
framework of our model, if the shock is located at x s = 100, t oscn ~ t a d v = 0.2s 
(vqpo = 5 for a M = 5M black hole. For a massive black hole, M = 10 8 M e , the 
time scale of QPO behavior could similarly vary from a few hours to a few days. A 
firm confirmation of this oscillation may be the test of accretion onto a very compact 
object through a shock wave. Exactly similar mechanism may be responsible for QPOs 
in disks around neutron stars as well. How long the oscillation will last depends on how 
sharp the resonance is. In a dwarf nova type outburst, accretion rate is swept through a 
large range, and therefore it is expected that in some range the resonance will occur. In 
white dwarfs candidates, such as SS Cygni, VW Hyi, and U Gem oscillations of 7 — lis 
(Mauche, Raymond k Mattei, 1995), ~ 14s (Van der Woerd et al, 1987) and 25 - 29s 
(Cordova et al. 1984) respectively have been seen. Another prediction of the model is 
that the time-period of oscillation (ss cooling time) should decrease (roughly linearly) 
as the accretion rate is increased. Similarly, oscillation may disappear completely if 
the accretion rate is out of the resonance width. These are clearly reported in recent 
observations (Mauche, 1995). Thus, we believe that the observational results may have 
agreed with our resonance model of oscillation. 

If the hard X-rays in galactic black hole candidates are actually produced in the 
post-shock region (Chakrabarti & Titarchuk, 1995), and the QPO oscillation is due to 
the oscillation of these shock waves, we see virtually no difference in the qualitative 
properties between the spectra of black hole candidates and that of neutron stars. 
The quantitative difference will come only from the presence (black hole) or absence 
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(neutron star) of the inner sonic point and thus from flow very close to the horizon. 
In other words, even if the flow becomes sub-Keplerian close to the black hole with 
very little radiation efficiency (Abramowicz et al. 1988, C89, C90, Chen & Taam, 1993, 
Chakrabarti & Titarchuk, 1995) it still finds a way, namely, through the shocks (C89, 
C90, CM95), to dissipate most of its energy at this 'boundary layer'. 

We have employed Smoothed Particle Hydrodynamics to illustrate this important 
time-dependent behavior of the flow. We observe that the numerical results agree with 
the analytical solution very accurately thus removing some of the misgivings attributed 
to SPH simulations. We have captured, quite successfully, very strong shock waves 
where Mach number jumps from M_ ~ 100 to M + ~ 0.45. We therefore believe that 
SPH could be a valuable tool to study astrophysical processes provided adequate care 
is taken to guarantee accuracy in interpolation method, which is accomplished here by 
the choice of very small particle size. 

We have studied the behavior of the accretion flow by considering only very 
simple physical processes such as compressional heating and bremsstrahlung cooling. 
The phenomenon that we discover, however, need not be confined to these systems 
alone. The oscillating behavior is due to generic physical processes such as heating and 
cooling and should be present even in complex situations where more than one heating 
and cooling mechanisms, such as Comptonization, line cooling etc. are present. We 
have shown in our simulation of radiative shocks in thick disks, that the oscillations are 
smoothed out to some degree and perhaps more closely behave like QPOs than periodic 
oscillations. In these simulations we show the existence of the hot corona which might 
be responsible for the hard component observed in galactic and extragalactic black hole 
candidate spectra. 

After completion of this work, we became aware of a very recent work of a purely 
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hydrodynamic simulation (Ryu et al, 1995) of a non-accreting, quasi-spherical flows 
around a Newtonian star. This simulation also shows some type of oscillation even 
when no cooling is added. These results seem to be due to the presence of a strong 
centrifugal barrier of a non-accreting high angular momentum flow. We have been able 
to reproduce these oscillations when both the cooling and the accretion are absent. In 
presence of accretion these oscillations go away and a steady state shock is achieved as 
in Molteni, Lanzafame k Chakrabarti (1994). Our present paper shows that to obtain 
a sustained periodic oscillation in an accreting disk, one requires to include cooling 
processes as in a realistic system. 

SKC would like to thank the hospitality of the Landessternwarte, and the Institute 
of Theoretical Astrophysics, University of Heidelberg where a part of the work was 
carried out. DM would like to thank the hospitality of the Institute of Theoretical 
Astrophysics, University of Heidelberg. Work of SKC at the Goddard Space Flight 
Center, is supported by a Senior Research Associateship award of the National Academy 
of Science. 

REFERENCES 

Abramowicz, M.A., Czerny, B., Lasota, J. P. k Szuskiewicz, E., ApJ, 332, 646 
Babul, A., Ostriker, J. P., k Meszaros, P. 1989, ApJ 347, 59 
Chakrabarti, S.K. 1989, ApJ, 347, 365 (C89) 

Chakrabarti, S.K. 1990, Theory of Transonic Astro-physical Flows (Singapore: World 
Scientific) (C90) 

Chakrabarti, S.K., k Molteni, D., ApJ, 1993, 417, 617 (CM93) 
Chakrabarti, S.K., k Molteni, D., MNRAS, 1995, 272, 80 (CM95) 
Chakrabarti, S.K., k Titarchuk, L.G. 1995, ApJ (in press) 
Chakrabarti, S.K., k Wiita, P.J., 1994, ApJ, 434, 518 
Chang K.M., k Ostriker, J. P., 1985, ApJ, 288, 428 

21 



Chen, X.M., k Taam, R.E., 1993, ApJ, 412, 254 
Chevalier, R.A. k Imamura, J.N., 1982, 261, 543 
Clavel, J. et al. 1990, MNRAS, 246, 668 
Clavel, J. et al. 1992, ApJ, 393, 113 
Cordova, F.A. et al. 1984, ApJ, 278, 739 

Dotani, T. et al. 1992, in Frontiers of X-ray Astronomy, (Universal Academy Press, 
Tokyo) p. 151 

Falle, S.A.G.E, 1995, in Shocks in Astrophysics, (Kluwar Academic Publications), (in 
press) 

Gingold, R.A. k Monaghan, J.J., MNRAS, 1977, 181, 375 

Imamura, J.N., Wolff, M.T. k Durisen, R.H., 1983, ApJ. 276, 667 

Kazanas, D. k Ellison, D.C., 1986, ApJ, 304, 178 

Langer, S.H. Chanmugam, G. k Shaviv, G., 1981, ApJ, 245, L23 

Mauche, C.W., Raymond, J.C. k Mattei, J. A., 1995, ApJ, (In press) 

Mauche, C.W. in Astrophysics in the Extreme Ultraviolet, ed. S. Bowyer k B. Haisch, 

(Camdridge: Cambridge University Press), in press 

Miller, J.S. k Peterson, B.M. ApJ 1990, 361, 98 

Monaghan J.J., Comp. Phys. Repts., 1985, 3, 71 

Monaghan J.J., Ann. Rev. Astron. Astrophys., 1992, 30, 543 

Molteni, D., Lanzafame, G., k Chakrabarti, S.K., ApJ, 425, 161 (MLC94) 

Molteni, D., k Sponholz, H. 1994, in Journal of Italian Astronomical Society, vol. 65-N. 

4-1994, Ed. G. Bodo k J.C. Miller 

Narayan, R., k Yi, I. 1994, ApJ, 428, L13 

Narayan, R., k Yi, I. 1995, ApJ (in press) 

Paczyhski, B, and Wiita, P.J. 1980, A&A, 88, 23 

Perola, G. C. et al. 1986, ApJ, 306, 508 

Peterson, B. M. et al. 1991, ApJ, 368, 119 



22 



Ryu, D, Brown, J., Ostriker, J. k Loeb, A, 1995, ApJ, in press 

Sponholz, H. k D. Molteni, 1994, MNRAS, 271, 233 

Shakura, N. I. k Sunyaev, R. A. 1973, A&A, 24, 337 

van der Woerd, H. k Heise, J. 1987, MNRAS, 225, 141 

Veilleux, S. k Zheng, W., 1991, ApJ. 377, 89 

Zheng, W., Veilleux, S. k Grandi, S.A., 1991, ApJ, 381, 418 



23 



Figure Captions 



Fig.l: Comparison of the time dependent shock locations as functions of time for a = 
0.4 (long-dashed line), a = 0.5 (solid line) and a = 0.75 (short-dashed line). Same 
accretion rate and angular momentum is used at the outer edge for all the three cases. 

Fig. 2: Comparison of the time dependent (long-dashed, short-dashed and dotted 
lines) and the steady state solutions (solid line) at various phases of oscillation when 
bremsstrahlung emission is considered. The time dependent solution oscillates close 
to the analytic solution. The dash-dotted curve represents the solution without any 
cooling effects. 

Fig. 3(a-b): The variations of (a) temperature (in units of 10 9 K) and (b) total bremsstrahlung 
loss (in arbitrary units) in various phases of oscillation. 

Figure 4: The loss of energy due to cooling as a function of time. Solid curve represents 
sinusoidal bremsstrahlung loss (a = 0.5) and the dashed curve represents the damped 
oscillation for a = 0.75 cooling. 

Fig. 5(a-d): Results of numerical simulations in three dimensional axisymmetric flows. 
In the upper half of each Figure, we have contours of constant Mach number M = 1 
(thick curves) and the velocity vectors and in the lower half we show the density contours 
(in units of density at the outer edge), (a) No cooling and (b-d) Post-shock corona at 
different phases of oscillations when the bremsstrahlung is included. 

Fig. 6: Oscillation of the net loss of energy due to bremsstrahlung as a function of time 
from a thick corona. 
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